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We derive the stationary probability distribution for a non-equilibrium system composed by an 
arbitrary number of degrees of freedom that are subject to Gaussian colored noise and a conservative 
potential. This is based on a multidimensional version of the Unified Colored Noise Approximation. 

By comparing theory with numerical simulations we demonstrate that the theoretical probability 
density quantitatively describes the accumulation of active particles around repulsive obstacles. In 
particular, for two particles with repulsive interactions, the probability of close contact decreases 
when one of the two particle is pinned. Moreover, in the case of isotropic confining potentials, the 
radial density profile shows a non trivial scaling with radius. Finally we show that the theory well 
approximates the “pressure” generated by the active particles allowing to derive an equation of state 
for a system of non-interacting colored noise-driven particles. 


Introduction 

A generic system, that is in thermal equilibrium at a 
temperature T, will be found in the neighbourhood of 
a configuration of energy E with a probability density 
given by the the Boltzmann factor exp[— E/ksT] [T) 2] - 
As an example, two Brownian colloidal particles, inter¬ 
acting through a conservative attractive force, will show 
an increased probability density for the low energy bound 
state. A quite different behaviour is observed in active 
particle systems [3j. Generally speaking, active matter 
is composed by biological or synthetic objects that are 
capable of absorbing energy from the environment and 
convert it into different kinds of persistent motions. Even 
when stationary states are reached, the probability dis¬ 
tributions can display large deviations from their equilib¬ 
rium Boltzmann counterparts. Those deviations are not 
just a matter of quantity, but a radically different qual¬ 
itative behaviour may be observed, like the widespread 
tendency to accumulate around repulsive objects. This 
“attraction for repulsion” is responsible for phenomena 
like, particle accumulation at solid walls [J] or the for¬ 
mation of bound states between repulsive objects U 15] - 
Our intuitive notion that particles like to stay where ex¬ 
ternal forces attract them to is biased by our familiarity 
with equilibrium statistical mechanics and needs to be 
replaced by novel statistical mechanics concepts that are 
capable to describe the stationary probability distribu¬ 
tions in systems of interacting active particles. In this 
context some schematic models have been proposed to 
model the dynamics of active particles as, for example, 
the “run and tumble” (RnT) model. The RnT dynam¬ 
ics is appropriate to describe the motion of bacteria such 
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as E. coli [THIS that swim along almost straight runs 
interrupted by random reorientations. In the case of col¬ 
loids propelled by chemical reactions the “active Brow¬ 
nian” model describes the motion of particles pushed 
by a force of constant magnitude that gradually reori¬ 
ent by rotational Brownian motion mHHi- However, 
despite the simplicity of the dynamics of these systems, 
it is hardly possible to find the analogue of the Boltz¬ 
mann distribution. Indeed, the Boltzmann prescription 
assigns a precise weight to a given configuration of po¬ 
sitions and momenta of particles at equilibrium. These 
particles are embedded in a space of dimensionality d , 
are subject to arbitrary external fields and mutually in¬ 
teract via whatsoever potential m • On the contrary, 
in the case of active particles the exact stationary prob¬ 
ability distribution is known only in rare instances as, 
for example, in the 1-dimensional RnT model in an ex¬ 
ternal force field [9]. The impossibility of writing explic¬ 
itly the stationary probability density prevents one from 
applying the standard methods of statistical mechanics. 
A Gaussian colored-noise model (GCN) can be used to 
reproduce the dynamics of passive colloidal particles im¬ 
mersed in a bath of dense swimming bacteria pa nsj. 
This model has been intensively studied in the past as 
the simplest model that could elucidate the basic physics 
of systems subject to time-correlated noise. Interestingly 
GCN was originally used to interpret the behaviour of 
very different physical systems such as noisy electronic 
circuits m and dye-laser radiation |T8l . The analytical 
study of GCN-driven systems resulted very challenging 
and led to the development of different approximation 
schemes aimed to reduce the complexity of the GCN- 
model to a tractable level pa cm. Among these ap¬ 
proaches one emerges by having a number of advantages 
with respect to the others. This is the Unified Colored 
Noise Approximation (UCNA) developed by Hanggi and 
Jung [20j that, under certain conditions, describes both 
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the small and the large correlation-time regimes, both in 
the high and low-friction limit [21] . More importantly for 
the present work the UCNA scheme can be generalized 
to a phase space of arbitrary dimensionality }22j . In this 
work we report, for the first time, the explicit formula of 
the stationary probability (obtained within the UCNA) 
for a system that is subject to a generic conservative po¬ 
tential and that is composed by an arbitrary number of 
degrees of freedom. We name this “multidimensional uni¬ 
fied colored noise approximated stationary probability” 
(MUCNASP). The MUCNASP plays basically the same 
role as the Boltzmann distribution for the approximated 
GCN-driven system. We show how the MUCNASP al¬ 
lows to predict several non-equilibrium properties of the 
active system in experimentally relevant cases where a 
simple external potential acts on a small number of de¬ 
grees of freedom. We focus on the case of steep repulsive 
interactions and spherically symmetric external poten¬ 
tials. In all these situations we use numerical simulation 
to test the quality of the approximation and find that the 
GCN-driven and the RnT particles display a strikingly 
similar behavior. In particular we show how our approxi¬ 
mated probability density captures very well the accumu¬ 
lation of the active particles around repulsive obstacles. 
Moreover the theory describes well the dependence on 
dimensionality of the probability density function when 
the active particles are confined by a circular repulsive 
wall. Understanding how the concept of pressure gener¬ 
alizes to active matter has become recently the subject 
of intense theoretical research [23lf25] . In this context we 
show how our theoretical probability density allows us to 
derive the pressure that the active particles exert on the 
repulsive walls and leads to the derivation of equations 
of state for the non-interacting active particle system. 
Finally we discuss the most relevant limitations of the 
present theory and suggest new routes to follow in the 
theoretical study of active matter. 


Results 

We consider the following set of stochastic differential 
equations: 


x = — V$ + Tj (1) 

where the position variables x = (aq,..., Xn) are deter¬ 
mined by the deterministic velocity — Vd* generated by 
a conservative potential < I>(x) and by a set of stochas¬ 
tic processes r/ = ( 771 , tjn). We assume that these 
are N independent Gaussian processes with zero mean 
(rij) = 0 and exponential time-correlation: {r]i(t) r lj{ s )} = 
— e _ l t_s l/ T . Here D is the diffusion coefficient charac- 
terizing the amplitude of the noise and r is its relaxation 
time. Note that here we absorb the mobility /i in the 


velocity field — V$ = /if, where f = —VC/ is the deter¬ 
ministic force generated by the potential energy function 
C/(x) = <h(x)//i. By using the UCNA Eq. © reduces to 
the (Stratonovich) Langevin equation [22] : 


[I + tH(x)]x = -V$(x) + D 1/2 T (2) 

where I = 6ij is the identity matrix, H(x) = d Xi d x ,$(x) 
is the Hessian associated with the potential <f>, and T is 
a set of independent white-noise sources having {Tj) = 0 
and (r i (t)r j -(s)) = 25ijS{t — s). We have found that, in 
the flow-free case, the steady state probability of finding 
the dynamical system of Eq. © in a specific configura¬ 
tion x is always proportional to the weight: 


f2(x) = exp 


$(x) 

D 


T 

2D 


|V$(x)| 2 


l|I + rH(x)|| 


(3) 


where |...| represents the norm of a vector and ||...|| in¬ 
dicates the absolute value of the determinant of a ma¬ 
trix. We have demonstrated the validity of Eq. © by 
deriving the corresponding Fokker-Planck equation from 
Eq. © [26] and solving it in the zero-current case by 
using the Jacobi’s formula [27) (see Supplemental Mate¬ 
rial) . 


One single degree of freedom 

As Eq. © is used for specific choices of the potential 
it reveals several interesting non-equilibrium properties 
of the active system under study. We initially focus on a 
simple one-dimensional case and study GCN-driven par¬ 
ticles when they are subject to a steep repulsive potential 
of the form $( 2 :) = Ax~ 12 setting A = 1. Such a poten¬ 
tial can be thought as a repulsive obstacle that perturbs 
the dynamics of the particles [28]. To verify the qual¬ 
ity of the UCNA, we integrate numerically the stochastic 
equation of motion Q in the presence of such a poten¬ 
tial. To this aim we have implemented a code for Euler 
integration of Eq. © , to be executed on GPU where the 
dynamics of many independent particles can be simu¬ 
lated in parallel [29]. We consider several different values 
of 0.1 < r < Is and 0.1 < D < 100/mi 2 /s, ranges that 
cover the typical persistence times and diffusivities of col¬ 
loids in bacterial baths, swimming bacteria such as E. 
coli [301131] and of chemically self-propelled Janus-type 
particles [T5] . The size of the simulation box is chosen 
to be L = 20 /tin (with periodic boundaries located at 
±L/2). Fig. 0 a) shows that in equilibrium (r = 0 ) the 
probability density decreases rapidly before the core of 
the repulsive potential is reached, whereas in the GCN 
the distribution peaks substantially in a region where the 
potential is very high before vanishing at the core. Note 
that the specific choice of the constant A = 1 defines 





3 


the size of the repulsive “wall” created by the external 
potential. The thickness of this impenetrable region is 
about 2 /jm and it depends very little on the values of D 
and t considered since the potential is steeply repulsive. 
In the GCN-driven system the exact probability distri¬ 
bution is unknown but it can be approximated by Eq. ([3| 
that reduces, for a single degree of freedom, to the known 
form pjj]: 


exactly for the RnT model, in particular the stationary 
probability distribution of the RnT model in presence of 
two hard walls is composed by two Dirac deltas plus a 
constant and the expression of (|'f> , |) has the same form of 
the one found in the UCNA (see Supplemental Material). 

Two interacting particles in one dimension. 


fl(x) = exp 


>I>(a;) 

D 


2D 


l* , (*)| S 


|l + r$"(x)| (4) 


where the prime indicates differentiation with respect to 
x. The approximate probability, obtained by normalizing 
Eq. Q P(x) = tl(x)/ f dxfl(x ), is plotted in Fig. [l](a) as 
a dashed line and it is found to reproduce well the numer¬ 
ical distribution at two well separated values of D and 
r. Knowing the probability we can also compute all the 
average quantities of interest, such as the average value 
of the modulus of the velocity (| < & , |) = f dxP(x) |$'(x)|. 
The theoretical (approximated) (|<fr , |) is compared with 
the numerical value in Fig. [ljb) where one sees that the 
(|<E> , |) UCNA prediction is very close to the numerically 
result at all values of D and r here investigated. The 
behaviour of the P(x ) and of (|$'|) can be qualitatively 
understood by considering the strong peaking of the n(x) 
close to the repulsive barrier. When the potential is 
very steep, as in the case of $ = x~ 12 , the maxima 
of Q(x) are found at x = x* where |$'(a;*)| ~ ^JD/t. 
It is clear that the probability peaks where the exter¬ 
nal potential balances the root mean-squared velocity 
of the particle induced by GCN (ij 2 ) = yjD/r. In 
this hard-wall limit is possible to approximate Q(x) in 
the neighborhood of x* by a strongly peaked function 
k(x) whose integral is « \JDt , as found by a saddle- 
point approximation, while far away from x* Eq. Q re¬ 
duces to unity (see Supplemental Material) and we can 
write: U(x) ~ k(x — x*) + k(x + x*) + 1. Note that 
\/Dt corresponds to the typical correlation length of the 
active motion. Integrating from —L/2 to L/2 we find 
the average velocity (|$'|) = 2 D/(L* + 2\/Dt), where 
L* = L — 2x* is the overall length available to the par¬ 
ticles and report this result as a dashed line in Fig. |TJ)b) 
where it captures the trend of (|<f>'|) obtained in simula¬ 
tions with the potential <f> = x~ 12 . By considering the 
force exerted by the particles located only on the right- 
hand side of the potential (x > 0) we find the force f x> o: 
(f x > o) = fi~ 1 D/(L* + 2 VDt), which corresponds to an 
equation of state being f x>0 the 1-dimensional pressure 
that the particles exert on the “wall” represented by the 
external potential. Note that if we consider N indepen¬ 
dent particles, in the limit r —► 0, and set D = /x ksT we 
arrive to the ideal gas law in Id: ( f x>0 } = NksT/L* and 
that this equilibrium value constitute an upper bound for 
the pressure of the active system (see dashed-dotted line 
in Fig. 0b)). Interestingly these results can be derived 


Up to this point we have derived results from the 
known formula Q considering only one degree of free¬ 
dom. We now use Eq. ([2]) to characterize the steady 
state properties of two interacting GCN-driven particles 
moving in Id. We consider the positions aq and x 2 of 
two particles interacting via a pair potential that is a 
function of the distance Ax = |aq — x 2 | between the par¬ 
ticles: $(xi,x 2 ) = 4>(Ax). We further assume that the 
particles are free to move in a Id space of extension L 
with periodic boundaries located at ±L/2. In this case 
Eq. © can be used to compute the probability of finding 
the two particles separated by Ax: 


II(Ax) = exp 


$(Ax) 

D 


2 T 

2D 


|$'(Ax)| 2 


|l + 2r$"(Ax)| 


(5) 

which is identical to Eq. Q with r replaced by 2r. 
This is at variance with equilibrium statistical mechan¬ 
ics in which the probability of finding two particles at 
a given distance does not vary if one of the two par¬ 
ticles is pinned at some fixed position [32HM1- To be 
more specific let us consider a repulsive potential of 
the form 4>(xi,x 2 ) = (xq — x 2 ) -12 . In tins case again 
the probability is well described by the MUCNASP as 
shown in Fig. 0b). Similarly the deterministic veloc¬ 
ity component experienced by one particle (| < I >, |) result¬ 
ing from simulations is well approximated by the theory 
(see Fig. [ljc) ). Again, in the limit of an infinitely steep 
potential (see Supplemental Material), we find that the 
area of the peak of Cl(Ax) is approximated by \J2Dt and 
(|<U|) = 2D/(L* + 2\/2Dt) which is plotted in Fig. [ljc) 
as a dashed line. This can be physically interpreted as 
follows: when both particles are free to move they can be 
more often found in contact since they move coherently 
in the same direction, this happens without the particles 
pushing onto each other, yielding a lower value of the 
average interaction force. It is important to note that 
these theoretical results cannot be derived analytically 
in the RnT model since the coupled dynamics of more 
particles makes the problem far too complicated. Nev¬ 
ertheless we have found that the MUCNASP produces 
results for the average (| < f>'|) that are very similar those 
found numerically for the RnT model despite the station¬ 
ary probability density has a very different form (see Sup¬ 
plemental Material). This suggests that the MUCNASP 
can be used as a convenient approximation also for calcu¬ 
lating the averages in RnT dynamics at least in the case 













4 


of steeply repulsive potentials. Note also that such a sce¬ 
nario is in agreement with the findings of Ref. 0 where 
it was demonstrated, by combining experiments and sim¬ 
ulations, that two colloids suspended in a bacterial bath 
tend to stay in contact because of the colored-noise forces 
induced by swimming bacteria. 






FIG. 1: (a) Probability density function of the position of a 
single GCN-driven particle in presence of the external po¬ 
tential $ = x~ 12 (shaded area). Full lines: simulations, 
dashed lines: theory, dashed-dotted line: Boltzmann dis¬ 
tribution. The curve with the higher peak corresponds to 
D = 100/rm 2 /s, r = Is and the one with the lower peak to 
D = 0.4/xm 2 /s, t = 0.1s (zoomed in the inset) (b) Average 
value of |$'| as a function of D for three different values of 
r = 0.1,0.325, 1 s from top to bottom respectively. Points: 
simulations, full lines: theory, dashed lines: theory in the 
limit of a hard potential, dashed dotted line: white noise case 
with a hard potential, (c) Probability density function of 
the distance between two GCN-driven particles interacting 
via the potential $ = A:r~ 12 , same legend as Fig. (a), (d) 
Average value of I®'! as a function of D for three different 
values of r = 0.1, 0.325,1 s from top to bottom respectively, 
same legend as as Fig. (b). 


Radially symmetric potentials. 

When the d-dimensional potential is spherically sym¬ 
metric, i.e. $(x) = < f>(r), Eq. ([3| simplifies to 


f l(r) = 0 exp 


<h(r) 
I) 

T$ , (r )| d_1 


2D 




|1 + T<f >,/ (r) 


( 6 ) 


where r 2 = J2i= i x i 2 an d 0 is the d-dimensional solid 
angle. Note that the Boltzmann distribution, obtained 
by setting r = 0 in Eq. ©, depends on the dimensional¬ 
ity only via the trivial term r d_1 while in the GCN-case 



FIG. 2: (a) Probability density function of the radial distance 
of one GCN-driven particle in presence of the spherically sym¬ 
metric external potential <I> = (r — R)~ 12 in 2d (shaded area). 
Full lines: simulations, dashed lines: theory, dashed-dotted 
line: Boltzmann distribution. The curve with the higher peak 
corresponds to D = 100/im 2 /s, t = Is and the one with the 
lower peak to D = 0.4 /xm 2 /s, r = 0.1 s (zoomed in the inset) 
(b) Average value of |<3? / 1 as a function of D for three different 
values of r = 0.1, 0.325,1 s from top to bottom respectively. 
Points: simulations, full lines: theory, dashed lines: theory in 
the limit of a hard potential, dashed dotted line: white noise 
case with a hard potential. 


this dependence is more complicated. To understand this 
issue, we consider GCN-driven particles in d = 2 in the 
presence of a circular repulsive potential of radius R of 
the form 4>(r) = (r — R)~ 12 where r = y/x 2 + y 2 and 
R = 5 pm. Simulation results show that particles ac¬ 
cumulate near the ring at r = R and the theoretical 
probability reproduces well this behaviour (see Fig.[2ja)). 
From Eq. © we can compute the averages of interest as 
the radial component of the velocity field (| < f>'|) and com¬ 
pare it with simulation results in Fig. [2jd) showing a good 
agreement. In the limit of an infinitely steep potential 
we have that f2(r) strongly peaks where |$'| ss y/D/r 
and reduces to unity elsewhere. The area of the peak 
can be approximated by 2itVDt\R* + \JDt | (see Sup¬ 
plemental Material), where R* = R — r* is the radial 
coordinate of the peak with r* ss 1 in the D-t range con¬ 
sidered. For the average radial velocity component we 
get (|$'|) = 2 D(R* + yfDr)/(R* 2 + 2 Dt + 2R*yfDr ) 
which is plotted as dashed lines in Fig. ©d) and fol¬ 
lows nicely the trend displayed by the numerical data. 
This is practically a colored-noise version of the ideal 
gas law in a circular container. This is clear if we set 
r = 0 to obtain (|4>'|) = 2 D/R* which is proportional to 
the average radial force (/) = 2 D/(R*p). Dividing this 
by the 2d “surface” 2ttR* we get the ideal gas pressure 
p = N(f)/(2nR*) = NksT/^nR* 2 ) for N independent 
particles and D = phsT. The unperturbed dynamics of 
the R.nT model in 2d is well understood ^S]i while the 
problem of a 2d symmetric potential is not tractable ana¬ 
lytically. However we have found that the simulation re¬ 
sults for the (|$ , |) in RnT model are very similar to those 
produced by the MUCNASP (see Supplemental Mate¬ 
rial) suggesting that MUCNASP could actually describe 
the behaviour of a wider class of active particles. An ex- 
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perimental situation, that is close to the spherically sym¬ 
metric case treated here, has been recently studied [36] . 
In these experiments it was observed a marked accumu¬ 
lation of swimming bacteria at the border of spherical 
liquid droplets. In this kind of experiment it would be in¬ 
teresting to check whether, in the dilute regime, the num¬ 
ber of bacteria found in contact with the border scales 
with the droplet radius and the characteristic run length 
as predicted by the MUCNASP (Eq. Q) in the spherical 
case. 


Discussion 

By using the unified colored noise approximation, we 
have derived the explicit formula for the non-equilibrium 
stationary probability (MUCNASP, Eq. (§) of a system 
composed by an arbitrary number of degrees of freedom 
subject to GCN. We have focused onto the case of steep 
repulsive potentials where the probability distribution of 
one single active particle tends to concentrate on the re¬ 
pulsive part of the potential oppositely to the case of a 
Brownian particle. Moreover we have verified that, as 
predicted by the MUCNASP, two active particles inter¬ 
acting repulsively behave differently with respect to equi¬ 
librium and are found in contact more often than if one 
particle is fixed. Finally the MUCNASP predicts that, 
when active particles are confined inside a repulsive ring- 
shaped barrier, the probability peaks on the boundary 
and the area of this peak increases with increasing ra¬ 
dius and with increasing persistence length. Surprisingly 
the results obtained by the MUCNASP are very close to 
those obtained for RnT particles for which an analytical 
solution is not available. As discussed before nnmum, 
the UCNA is accurate in those portions of phase space 
where the all the eigenvalues of the Hessian matrix are 
positive. This restriction defines where Eq. (|3| can be 
used as a valid approximation for the probability of GCN- 
driven system. Moreover it appears difficult to derive an 
explicitly formula for probability including also Brownian 
fluctuations. Nevertheless, to our knowledge, the MUC¬ 
NASP is the only available explicit probability formula 
accounting for multiple active degrees of freedom and, 
provided that all eigenvalues are positive, becomes exact 
both in the limit of r —» 0 and r —> oo. This makes 
the MUCNASP a valuable schematic model for tackling 
the many-body problem in active matter. For example 
it would be interesting to check whether, at the mean- 
field level, the MUCNASP can predict a motility-induced 
phase separation mmm or the colored-noise induced 
shift in the synchronization threshold of the noisy Ku- 
ramoto model [55] , 
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